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Abstract. We present a program tool, SimClust, designed for Monte-Carlo 
modeling of star clusters. It populates the available stellar isochrones with stars 
according to the initial mass function and distributes stars randomly following 
the analytical surface number density profile. The tool is aimed at simulating 
realistic images of extragalactic star clusters and can be used to: (i) optimize 
object detection algorithms, (ii) perform artificial cluster tests for the analysis 
of star cluster surveys, and (iii) assess the stochastic effects introduced into 
photometric and structural parameters of clusters due to random distribution 
of luminous stars and non-uniform interstellar extinction. By applying Sim- 
Clust, we have demonstrated a significant influence of stochastic effects on 
the determined photometric and structural parameters of low-mass star clus- 
ters in the M 31 galaxy disk. The source code and examples are available at the 
SimClust website: http://www.astro.ff .vu.lt/software/simclust/. 

Key words: galaxies: star clusters - galaxies: individual (M31) - methods: 
numerical - techniques: photometric 

1. INTRODUCTION 

Recent studies of star clusters have extended our understanding of these grav- 
itationally bound systems, which have also been successfully used to investigate 
evolution of their host galaxies; see a review by Kroupa (2008) and references 
therein. The recent observational data (see review by Renzini 2008) suggest the 
presence of multiple stellar populations in star clusters. Theoretical modeling 
(e.g., Lamers et al. 2006; Kruijssen & Lamers 2008) indicates that the stellar 
mass function evolves from the initial mass function (IMF) due to dynamical ef- 
fects as the cluster ages. However, the assumption that a star cluster consists of a 
simple stellar population (SSP) is still an attractive simplification in extragalactic 
studies. 

Although due to the improved quality of observations the detection of extensive 
low-mass extragalactic star cluster samples is now possible, the straightforward ap- 
plication of SSP models to derive their evolutionary parameters is valid only for 
high-mass clusters (Cervirio & Luridiana 2004), even neglecting effects of dynam- 
ical evolution and loss of stars. The morphology of young low-mass star clusters 
also depends on the stochastic distribution of luminous stars and non-uniform ex- 
tinction, therefore, their structural, as well as photometric parameters are biased. 
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These stochastic uncertainties are also very important for the study of selection 
effects in cxtragalactic surveys of unresolved and semi-resolved star clusters, which 
should be quantified in terms of age, mass, extinction and size. 

Therefore, we have developed a program tool, SimClust, by adopting Monte- 
Carlo modeling, which was previously used to study cluster color distributions 
(e.g., Girardi & Bica 1993; Bruzual 2002), but extending this approach to simu- 
late realistic images. Thus, proper investigation of stochastic effects and observa- 
tional errors on structural parameters and aperture CCD photometry, as well as 
the effects of non-uniform extinction in a star cluster, become now possible. The 
program code is available under the GNU licenseQ and may be modified to incor- 
porate different isochrones and stellar surface number density distribution profiles. 
The source code and examples are available at the SimClust websitfl 

SimClust was applied to investigate stochastic effects and the validity of the 
SSP model fitting for compact star cluster sample in the disk of the M 31 galaxy 
(Narbutis et al. 2008). In this study we have simulated numerous images of 
artificial clusters of typical age, mass and size, characteristic of the M 31 clusters 
(Vansevicius et al. 2009). Significant influence of stochastic effects was found for 
low-mass (<;3000 Mq) star clusters younger than ~ 100 Myr. 

This program tool shares a similar approach to the star cluster simulation prob- 
lem investigated by Popescu & Hanson (2008). SimClust uses different stellar 
isochrone sets, initial mass function sampling technique and modeling of inter- 
stellar extinction, thus the user is provided with an independent tool to perform 
Monte-Carlo simulations of extragalactic star clusters and to check the robustness 
of cluster models. 

We present the star cluster simulation method in § 2, describe the SimClust 
program in § 3 and investigate stochastic effects for cluster sample in the M 31 
galaxy in § 4. We also discuss the implications of Monte-Carlo modeling in extra- 
galactic star cluster studies on both the estimation of star cluster sample complete- 
ness and determination of object parameters in the stochastic model framework. 

2. STAR CLUSTER SIMULATION 

SimClust program populates a stellar isochrone of a given age and metallicity 
with stars. Randomly generated masses of stars are weighted on the IMF and the 
process is repeated until a given mass of model star cluster, M, is reached. Stellar 
magnitudes are computed assuming a distance to the artificial "extragalactic" 
cluster, and stars are scattered in the image plane by the analytical surface number 
density profile. The effect of non-uniform extinction can be applied for individual 
stars. Finally, realistic multi-band images of model cluster are rendered. We 
describe the details of each step below. 

2.1. Isochrones and extinction - stellar magnitudes 

For stellar magnitudes in different photometric passbands as a function of mi- 
tiafi^ stellar mass in solar-mass unit, m, age, t, and metallicity, [M/H], we use 
a new set of PadovE@ isochrones by Marigo et al. (2008). These isochrones are 

1 GNU General Public License. 

2 Download SimClust from http://www.astro.ff .vu.lt/software/simclust71. 

3 Initial and actual (smaller due to mass loss) stellar masses are provided in isochrone tables. 

4 Padova isochrones: http://stev.oapd.inaf.it/cgi-bin/cmd. 
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combined of models: (i) up to early- AGB from Girardi et al. (2000), (ii) detailed 
TP-AGB for m < 7M© from Marigo & Girardi (2007) and (iii) for m > 7M 
from Bertelli et al. (1994), and take into account many critical aspects of AGB 
evolution, important for star cluster simulations. In this study we have used the 
Marigo et al. (2008) isochrones in the Johnson-Cousins UB VRIJHK photometric 
system (Mai'z Apellaniz 2006; Bessell 1990). 

The isochrone consists of points for a discrete stellar mass with tabulated 
corresponding magnitudes for each passband. For every randomly generated star 
(see §2.2), we calculate its magnitude in each passband by linear interpolation 
between the closest mass points on the isochrone. The accuracy of the magnitude is 
directly related to the mass sampling on the isochrone, and the linear interpolation 
should have negligible effect on the accuracy of the integrated cluster magnitude. 

The interstellar extinction is applied for individual stars assuming the standard 
extinction law (Cardelli et al. 1989), calculated as a function of the given color 
excess, Eb-v- To investigate the effects of non-uniform (differential) extinction in 
young clusters (Yadav & Sagar 2001), we assume for simplicity that reddening 
of individual stars in the cluster is described by a Gaussian distribution over 
mean Eb~v and a standard deviation o~e b -v This distribution is generated 
using gsl_ran_gaussian function from the GNU Scientific Librarjo- Setting the 
°~e b -v = 0j all stars are reddened by the same amount. In case of oe b -v > 0, the 
lower limit of reddening, corresponding to the foreground extinction, is applied. 

2.2. IMF - stellar mass distribution 

The stellar initial mass function defines a number of stars, dN, within a stellar 
mass range, dm. In this work we use the universal IMF (Kroupa 2001), defined as 
the multi-power law, which has the general form of 

f (m) = dN/dm = Cimf ■ h ■ m"> , (1) 

where cti is the IMF slope for the mass interval [rai_i,mj, Cimf and bi are the 
normalization and function continuation constants, respectively. 

We follow the practical numerical formulation of the IMF by Pflamm-Altenburg 
& Kroupa (2006), and implement their shared library "Libimf'0 in our code; 
namely the function imf _init_multi_power, which allows us to define an arbitrary 
number i of mass intervals [m,_i,mj], and corresponding power law slopes cti. 

Because of the lowest stellar mass limit m m j n = 0.15 M Q for the isochrones 
used, the IMF slope (Eq. 2, Kroupa 2001) transforms into a simple equation: 

f a = -1.3, m/M e < 0.50; , . 

{ ai = -2.3, m/M > 0.50. { > 

Before a random number generator can be applied to sample the IMF, it has 
to be normalized, i.e., the constant Cimf has to be defined. We follow the normal- 
ization strategy by Weidner & Kroupa (2004) , implemented by Pflamm-Altenburg 
& Kroupa (2006) in the "Libimf" library as a function imf _norm_wk04. It requires 
two additional free parameters: (i) m max * = 150 Mq - the maximum physically 
possible stellar mass, and (ii) m max - the expected maximum stellar mass in a 

5 GSL version 1.9: http://www.gnu.org/software/gsl/ . 

6 Libimf: http://www.astro.uni-bonn.de/~webaiub/dwld/llbimf-05-10-06.tar.gz . 
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Fig. 1. Properties of a model star cluster with t = 100 Myr and "luminous" mass 
A/ium = 3000 Mq, consisting of JVl um = 5300 stars. The cumulative mass distribution 
function E(m) of the IMF (see Eq. [3| and the cumulative radial distribution function 
^(r) of the King model (see Eq. [6}, displayed in panels (a) and (b), respectively, are used 
to sample mass, m, and radial distance, r, by means of uniform random numbers, rnd. 
The stellar mass distribution, £(log(m)) = dN/d\og(m), is displayed vs. mass, m, in 
panel (c), overplotted with slope cto — —1.3 and ot\ — —2.3 lines of IMF (see Eq.0. The 
lower and upper mass limits are defined by the isochrone. The surface number density of 
stars, ip(r) = dN/dS, is plotted vs. radial distance, r, in panel (d), obeying King (1962) 
model profile of r c = 0.75 pc and r t = 15 pc (see Eq. [5]); half-light radius rh ~ 1.7 pc is 
indicated by an arrow. 



given cluster of mass M. A reasonable assumption is to define both Cimf and 
Ti max in such a way, that there is only one massive star of mass m max in the mass 
range [m max , rftmax*]. Thus, the expected maximum stellar mass is a function of 
cluster's mass, i.e., TO max = w max (M). The normalization is a required procedure 
for the present application of the Monte-Carlo method, since now the IMF can be 
treated as a probability density function in the range between m and m + dm. 

The IMF sampling technique can be summarized as follows: from Eq. [1] a 
cumulative mass distribution function, giving a probability that mass m? < to, is 
constructed: 

j a 

H(m) = / ^(m')dm'. (3) 



This function (Eq. [3|) is continuous between and 1 (see Figure la). Thus, a 
uniform random number, rnd, obtained from long periodicity gsl_rng_taus ran- 
domizer function of the GNU Scientific Library, is attributed to H(m) = rnd. The 
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stellar mass m values are sampled from the IMF via inverse function 

m = F^ 1 (5(m) = rnd) , (4) 

which has an analytical solution. We use the "Libimf's" imf _dice_star_cl func- 
tion to generate an array of stellar masses according to the supplied IMF. 

The mass of the cluster derived by the SSP model fitting and used in pop- 
ulation analysis is usually the initial mass at the cluster's birth. Following this 
definition, we generate the initial mass function of stars at age t = Myr, where 
individual stellar masses populate the range [m m i n , m max ], as defined previously. 
We denote the sum of initial stellar masses to the mass of the cluster, M. When 
the isochrone of age t > Myr is sampled, only stars, which are present at that 
age (m max = m max (t) < m max (t = OMyr)) contribute to the luminosity of cluster, 
i.e., "luminous" mass, M\ um , is always lower than the "initial" mass, M, of the 
cluster. Also, the mass-loss of stars due to stellar winds can be accounted for and 
directed to the output. 

An example of the resulting stellar mass distribution, £ (log (to)) = dN/ d log (to), 
displayed vs. mass, m, of t = 100 Myr model star cluster of "luminous" mass 
Mi um — 3000 M Q , consisting of N\ um — 5300 stars, is shown in Figure lc, over- 
plotted with the IMF slope a = —1.3 and ol\ = —2.3 lines (see Eq. [5]). 

2.3. Radial surface number density profile - 2-D spatial distribution 

Once the mass and multi-band magnitudes of a star are known, its spatial 
position has to be defined. For simplicity we assume that the "observed" positions 
of the stars are distributed randomly in the 2-D plane, obeying a circular King 
(1962) model. 

The King model describes the radial surface number density profile, i.e., the 
number of stars, dN, per unit area, dS and is defined by the central surface number 
density, tpo, the core radius, r c , and the tidal radius, rv 



ip (r) = dN/dS = dN/2-nr dr = xp 



-1/2 



-1/2 



(5) 



where r is the distance from cluster's center. 

The cumulative radial distribution function of the circular King model giving 
probability that the distance r' < r, is defined by formula: 



tt(r) =2^C Ki ng / r'^(r')dr' 



In (1 + (r/r c ) 2 ) - 4 y/(r*+r*)/\ - ^rJ/X) + r 2 / X 



(6) 



where A = r 2 + rf, and CKine is the normalization constant, defined as 



Cki: 



2tt J r'ip(r')dr' 



(7) 
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The vff (r) function (Eq. E|) is continuous between and 1 (see Figure lb). 
Thus, a uniform random number, rnd, is attributed to \& (r) = rnd to sample the 
distance r values from the King model via inverse function 



which docs not have an analytical expression. Since the function \&(r) is monotonously 
increasing, therefore, for a generated random number, rnd, we can solve Eq. [5]via 
Eq. El iteratively by the bisection method until the solution of r converges. 

An example of the resulting surface number density of stars, tp(r) = dN/dS, 
displayed vs. radial distance, r, of the model star cluster, consisting of N\ um = 5300 
stars, is shown in Figure Id, overplotted with the King (1962) model profile of 
?'c = 0.75 pc and ft = 15 pc (see Eq. [5]); half-lighi0 radius ~ 1.7 pc is indicated 
by an arrow. 

Finally, the azimuth angle, ip, is assigned randomly for a star. Stellar positions 
in parsecs (x, y) are computed in respect to the cluster center: 



2.4- Realistic images 

The distance to the star cluster, D, and the photometric zero-point are used 
to transform the absolute magnitudes of stars to the instrumental magnitude sys- 
tem in each photometric passband. Stellar positions in the image are computed 
assuming the provided pixel scale, 0, and the center coordinates of star cluster. 

Finally, these data are passed to the external "SkyMaker'H program, which 
renders realistic images for each photometric passband, assuming the provided 
point-spread function (PSF). Usually the PSF is derived from the real image by 
using the DAOPHOT program package (Stetson 1987) implemented in the IRAF 
program system (Tody 1993). Since "SkyMaker" has many image rendering op- 
tions, they can be used to create realistic image noise properties, or even simulate 
images for an arbitrary telescope configuration or observational conditions. 

3. SimClust - PROGRAM DESCRIPTION 

Here we present a short description of SimClust usage, the properties of input 
files, parameter configuration and output files. The schematic representation of 
the data flow is shown by means of a block diagram given in Figure 2. 

SimClust is written in the Perl and C++ languages and can be compiled 
under the "UNIX/Linux" systems. Installation script compile checks that exter- 
nal program "SkyMaker" and all necessary libraries are available on the computer 
and builds executables. The configuration file example, which contains main pa- 
rameters of the model star cluster and instructions on how the images should be 
simulated, is passed to SimClust which controls two sub-programs, ClsCreator 
and FitsCreator. 

7 In this case this is the radius containing the half number of the stars. 

8 SkyMaker version 3.1.0: http://terapix.iap.fr/soft/skymaker. 
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(8) 
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Fig. 2. Schematic representation of SimClust, which consists of two sub-programs 
ClsCreator and FitsCreator. Arrows indicate the parameter and data flow. 



The ClsCreator, which randomly populates the input stellar isochrones, has 
to be executed first. There are three possible output options: (i) .cmd file, con- 
taining magnitudes of individual stars, (ii) . sum - integrated magnitudes of star 
clusters, and (iii) .dat - stellar magnitudes and positions. Later the FitsCre- 
ator can be asked to pass .dat files and PSF to the "SkyMaker" to simulate 
cluster images (.fits). 

The detailed description of each configuration parameter is documented in the 
help file and illustrated with examples. 



4. STOCHASTIC EFFECTS OF STAR CLUSTERS 

The accuracy of cluster age and mass from broad-band photometry, derived 
by applying traditional SSP models, as well as the accuracy of structural cluster 
parameters, depend on the clusters age, mass and size, image spatial resolution, 
data analysis techniques and stochastic effects. To account properly for stochastic 
effects, we have developed SimClust and applied it to investigate properties of star 
clusters in the extensively studied South- West field of the M31 galaxy (Kodaira 
et al. 2004, 2008; Narbutis et al. 2006, 2007b, 2008). Previously, for this cluster 
sample the influence of individual semi-resolved bright stars on the accuracy of the 
structural cluster parameters was discussed by Sableviciute et al. (2006, 2007), 
and the influence on the accuracy of aperture photometry was studied by Narbutis 
et al. (2007a). 

In the following we provide an example of SimClust application to study star 
clusters in M 31. We describe the cluster simulations and the expected photometric 
and half-light radii uncertainties, derived for typical star clusters of our sample. 
For the present study we assume the typical parameter values of compact star 
clusters (Vansevicius et al. 2009): age ~ 100 Myr, mass ~3000M Q and half-light 
radius ~ 1.5 pc. 
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Fig. 3. Modeling of stochastic effects. Panel (a). The diagram V vs. B—V of 238 high 
probability star cluster candidates in M31 (gray circles) from Vansevicius et al. (2009), 
overplotted with five cases of model clusters for the ages t = 20, 100, 500 Myr, 2.5 
and 12.5 Gyr, and "luminous" masses log (Mium/iW©) = 3.5, 3.5, 3.5, 4.0 and 5.0 (blue, 
green, magenta, cyan and red symbols), 1000 models for each case. Asterisks indicate the 
position of the SSP model, vector indicates interstellar reddening of Eb-v = 0.25. Panel 
(b). The diagram U-B vs. R-I, symbols are the same as in (a). Panel (c). Distributions 
of model clusters half-light radii, measured using a curve-of-growth method, for the ages 
t = 20, 100 and 500 Myr; the dashed line indicates the true = 2.6 pc. 



4-1. Star cluster model setup 

Since the majority of objects in the M31 sample are young low-mass star 
clusters, we were motivated to investigate the influence of stochastic distribution 
of bright stars on both photometric and structural properties. For this purpose, 
we have chosen parameters adequate to the characteristic values of our sample: 
(i) five age cases, t = 20, 100, 500 Myr, 2.5 and 12.5 Gyr; (ii) "luminous" masses, 
which correspond to the adopted ages, log(Mi um /M©) = 3.5, 3.5, 3.5, 4.0 and 5.0; 
(iii) metallicity, [M/H] = -0.4 dex; (iv) reddening, Eb-v = 0.25; (v) distribution 
of stars according to the circular King model of r c = 0.75 pc and r t = 15 pc, i.e., 
rh ~ 1.7 pc. A thousand model clusters were produced for each age-mass cases. 

Model star clusters were placed at the distance of M31, D = 785 kpc, derived 
by McConnachie et al. (2005); thus 1" = 3.8 pc. For simplicity, zero photometric 
background and noise was assumed for image simulation. 
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4-. 2. Photometric properties 

To examine the photometric properties of model clusters, images equivalent to 
the Local Group Galaxies Survey data (Massey et al. 2006) with homogenized 
PSF of FWHM w 1.5" were simulated; image scale 6 = 0.27" per pixel. The 
UB VRI photometry was performed with IRAF's phot through apertures of 3" in 
diameter - a typical size used for the sample clusters (see Narbutis et al. 2008 for 
details of photometry). 

The resulting diagrams V vs. B-V and U—B vs. R—I are displayed in Fig- 
ures 3a and 3b, respectively, overplotted on 238 high probability cluster candidates, 
selected by Vansevicius et al. (2009). The stochastic scatter of model luminos- 
ity and color dominates over the typical photometric errors. The indicated SSP 
model points (Marigo et al. 2008) of [M/H] = -0.4 dex, reddened by E B ^ V = 
0.25, show that SSP colors at t — 20 Myr and 500 Myr do not represent model 
cluster colors and are disturbed asymmetrically around the SSP position in two- 
color diagram due to stochastic existence of bright blue and red supergiants, RGB 
and AGB stars, prominent at these ages. However, for the remaining three model 
ages, the centers of color distribution coincide well with the SSP models. The color 
distribution of sample clusters is covered by the extent of stochastic models with 
assumed typical reddening of Eb-v = 0.25. We note, that the scatter of colors 
of t = 100 Myr and log (Mi um /M Q ) = 3.5 models resemble that of real objects 
located around [0.4,-0.2] in Figure 3b. 

4-3. Structural properties 

To examine the structural properties of model clusters, images equivalent to 
the Suprime-Cam U-band data of FWHMpgp ~ 0.7" were simulated; the image 
scale 9 = 0.20" per pixel (see Narbutis et al. 2008 for Suprime-Cam survey data 
details). 

IRAF's phot was applied to obtain the profile of the enclosed flux of a cluster 
model using the curve-of-growth method, yielding the half-light radius, rZ. Note, 
however, that observational effects (PSF) were not taken into account. 

The derived half-light radii r^ data histograms of 10 3 models for ages t — 
20, 100 and 500 Myr were constructed. After smoothing, they are displayed in 
Figure 3c. The scatter around the true r h ~ 2.6 pc value is due to a stochastic 
distribution of bright stars in model star clusters. The stochastic effect is most 
significant for the t — 20 Myr models, and decreases to older and more massive 
clusters. The values of rj are smaller for clusters, which have bright stars residing 
in the cluster centers, dominating the overall luminosity. Contrary, when bright 
stars are located far from the center, the derived rjj values are increased. The 
standard deviation of r£ distribution of the youngest model clusters (t = 20 Myr) 
is ~ 0.5 pc and decreases to ~ 0.1 pc for the oldest ones (t = 500 Myr). 

5. DISCUSSION AND CONCLUSIONS 

Simulations of star cluster images are becoming important to evaluate the 
performance of analysis techniques of observational data, as have been shown by, 
e.g., Ascenso et al. (2008) for the case of mass segregation in star clusters. Popescu 
& Hanson (2008) emphasize the importance of realistic star cluster simulation with 



9 Padova SSP models: http://stev.oapd.inaf.it/cgi-bin/cmd. 



360 



V. Deveikis, D. Narbutis, R. Stonkute et al. 



incorporation of interstellar extinction to answer the question about the lack or 
existence of massive star clusters in galaxies. The current knowledge about the 
most massive clusters can be limited due to unknown selection effects. 

There is no technique available (Cerviho & Luridiana 2006) to properly account 
for stochastic effects when photometric data are interpreted in the framework of 
SSP models, except of star cluster simulation and assessment of the reliability of 
the SSP model performance in certain age and mass domains, ft might turn-out 
that the interpretation of properties of the cluster population, especially in the 
young low-mass domain, will change dramatically once the stochastic effects are 
taken into account. On the other hand, the fundamental limitation on the amount 
of information, which can be extracted from integrated photometry alone, will be 
established and other techniques shall be developed. 

It should be emphasized here, that even for the high-mass star clusters, where 
stochastic uncertainties are comparable to the observation photometric errors, pa- 
rameter degeneracies in broad-band photometry are unavoidable; see Bridzius et 
al. (2008) for the performance analysis of the UB VRIJHK photometric system. 

The tests performed in this study show that the stochastic effects are significant 
for low-mass young clusters. The age degeneracy due to stochastic color scatter 
is observed in all model cases. The reddening vector in Figure 3 implies that 
stochastic color shift could also mimic reddening. Therefore, when SSP model 
fitting is applied, the derived cluster mass becomes more uncertain, not only due 
to the miss-match of age and extinction, but also due to the intrinsic scatter 
of luminosity. Visual inspection of the model cluster color images shows a close 
resemblance to the real compact cluster candidates in respect to their semi-resolved 
appearance and stochastic emergence of bright blue or red stars. 

The morphology of low-mass star clusters is strongly dependent on the stochas- 
tic distribution of luminous stars and the spatial resolution of the image. A usual 
"by-eye" selection technique leads to systematic biases in the magnitude limited 
samples, thus automatic detection should be used instead. To properly estimate 
the completeness of cluster sample, the artificial star cluster tests should be per- 
formed using realistic simulations of cluster images on individual star basis, instead 
of just assuming the smooth surface brightness distribution profiles. 

The stochastic SSP model framework, which could replace the traditional SSP 
model fitting to quantify ages, masses and extinctions of star clusters from their in- 
tegrated broad-band photometry, is still to be developed. The SimClust program 
tool may serve as a basis for this method, as well as for extensive extragalactic 
star cluster studies. 
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